Raport z projektu zaliczeniowego

Opis danych i cel projektu

Zebrane dane zawieraja zbior informacji na temat pomiaru cech fizycznych krabow hodowlanych.

Zbior ten ma 9 roznych cech:

• Plec (M - Mezczyzna, K - Kobieta, I - Nieokreslony),

• Dlugosc kraba (w stopach),

• Srednica kraba (w stopach),

• Wysokosc kraba (w stopach),

• Waga kraba (w uncjach),

• Waga kraba bez skorupy (w uncjach),

• Waga samych wnetrznosci kraba (w uncjach),

• Waga skorupy kraba (w uncjach),

• Wiek kraba (w miesiacach).

Zestaw tych danych posiada 3893 obserwacje. Kazdy krab nalezy z zalozenia do jednej z trzech klas (podzial na dana plec). Zatem badamy 8 roznych cech.

Dane pozyskane zostaly z https://www.kaggle.com/datasets/sidhus/crab-age-prediction

Celem projektu jest przeprowadzenie klasyfikacji danych bez nadzoru i pod nadzorem na podstawie wyselekcjnowanych 9 parametrów.

Plan analizy:

1. Prezentacja statystyk sumarycznych, histogramow i wykresow gestosci

2. Analiza jednowymiarowa obserwacji odstajacych oraz wielomodalnosci

3. Analiza korelacji i redukcja wymiaru metoda PCA i UMAP

4. Analiza skupien - klastrowanie

5. Klasyfikacja pod nadzorem z uzyciem lasow losowych

6. Analiza z uzyciem LDA i QDA

1. Statystyki sumaryczne i histogramy

Prezentacja naszych danych:

head(krab)

##   Sex Length Diameter Height    Weight Shucked.Weight Viscera.Weight
## 1   F 1.4375   1.1750 0.4125 24.635715      12.332033       5.584852
## 2   M 0.8875   0.6500 0.2125  5.400580       2.296310       1.374951
## 3   I 1.0375   0.7750 0.2500  7.952035       3.231843       1.601747
## 4   F 1.1750   0.8875 0.2500 13.480187       4.748541       2.282135
## 5   I 0.8875   0.6625 0.2125  6.903103       3.458639       1.488349
## 6   F 1.5500   1.1625 0.3500 28.661344      13.579410       6.761356
##   Shell.Weight Age
## 1     6.747181   9
## 2     1.559222   6
## 3     2.764076   6
## 4     5.244657  10
## 5     1.700970   6
## 6     7.229122   8

cat('podsumowanie statystyczne','\n')

## podsumowanie statystyczne

apply(krab[,-1],2,summary)

##           Length Diameter    Height    Weight Shucked.Weight Viscera.Weight
## Min.    0.187500 0.137500 0.0000000  0.056699      0.0283495     0.01417475
## 1st Qu. 1.125000 0.875000 0.2875000 12.672227      5.3438808     2.66485300
## Median  1.362500 1.062500 0.3625000 22.792998      9.5396068     4.86193925
## Mean    1.311306 1.020893 0.3493739 23.567275     10.2073420     5.13654636
## 3rd Qu. 1.537500 1.200000 0.4125000 32.786197     14.2739732     7.20077300
## Max.    2.037500 1.625000 2.8250000 80.101512     42.1840560    21.54562000
##         Shell.Weight       Age
## Min.      0.04252425  1.000000
## 1st Qu.   3.71378450  8.000000
## Median    6.66213250 10.000000
## Mean      6.79584412  9.954791
## 3rd Qu.   9.35533500 11.000000
## Max.     28.49124750 29.000000

cat('odchylenie standardowe','\n')

## odchylenie standardowe

apply(krab[,-1],2,sd)

##         Length       Diameter         Height         Weight Shucked.Weight 
##      0.3004306      0.2482327      0.1049762     13.8912007      6.2752747 
## Viscera.Weight   Shell.Weight            Age 
##      3.1041331      3.9433917      3.2209666

cat('skosnosc','\n')

## skosnosc

apply(krab[,-1],2,skewness)

##         Length       Diameter         Height         Weight Shucked.Weight 
##     -0.6508601     -0.6171649      3.3130681      0.5187028      0.7028435 
## Viscera.Weight   Shell.Weight            Age 
##      0.5780672      0.6074819      1.1042846

Wniosek: Zmienna Weight ma zdecydowanie najwieksza wariancje oraz najwieksza wartosc srednia sposrod innych zmiennych. Oznacza to, ze ta zmienna ma duza zmiennosc danych. Natomiast najmniejsza zmiennoscia danych w badanej probce wykazuje sie zmienna Height, jednak widzimy w podsumowaniu statystycznym, ze jej wartosc maksymalna znacznie odstaje od pozostalych jej wartosci. Zmienne Length i Diameter maja ujemna skosnosc, co wskazuje na lewoskosnosc, natomiast pozostale zmienne sa prawoskosne.

for(i in 2:9) hist(krab[,i],main=colnames(krab)[i])

Wniosek: Z uzyskanych histogramow mozna wyczytac, ze zmienna Height ma na pewno odznaczajace sie wartosci odstajace, poniewaz jej histogram jest skupiony tylko na jednej stronie rozkladu. Wiadomo juz, ze ta zmienna ma duza skosnosc, co takze wplywa na asymetrie rozkladu. Zmienne Length i Diameter to rozklady lewoskosne, natomiast zmienna Age jest prawoskosna. Nie ma podejrzen co do rozkladu normalnego, jednak trzeba zbadac zmienne Height, Weight, Shucked.Weight, Viscera.Weight i Shell.Weight na wielomodalnosc.

for(i in 2:9) plot(density(krab[,i]),main=colnames(krab)[i])

Wniosek: Z badania gestosci nie wynika jasno, ktora zmienna sposrod wszystkich ma wiekszy potencjal jako zmienna rozrozniajaca pomiedzy plciami kraba, poniewaz wykresy te nie wskazuja ani na wielomodalnosc ani na rozklad normalny.

Mozemy zauwazyc pewna zaleznosc, mianowicie na histogramie oraz na wykresie gestosci zmienne Length oraz Diameter sa do siebie bardzo podobne. Zastosujmy wykres kwantyl-kwantyl.

qqplot(krab$Length,krab$Diameter,main='Wykres kwantyl-kwantyl')
abline(0,1)

Wniosek: Obie zmienne maja mocno podobny do siebie rozklad, poniewaz punkty ulozyly sie wzdluz linii referencyjnej. Mowi nam to, ze miedzy nimi jest duza korelacja. Jest to logiczne, poniewaz dlugosc i srednica kraba powinny byc scisle powiazane miedzy soba.

2. Analiza jednowymiarowa obserwacji odstajacych oraz wielomodalnosci

apply(krab[,-1],2,function(x) shapiro.test(x))

## $Length
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.96863, p-value < 2.2e-16
## 
## 
## $Diameter
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.96961, p-value < 2.2e-16
## 
## 
## $Height
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.88202, p-value < 2.2e-16
## 
## 
## $Weight
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.97309, p-value < 2.2e-16
## 
## 
## $Shucked.Weight
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.96325, p-value < 2.2e-16
## 
## 
## $Viscera.Weight
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.96909, p-value < 2.2e-16
## 
## 
## $Shell.Weight
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.97133, p-value < 2.2e-16
## 
## 
## $Age
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.93235, p-value < 2.2e-16

Wniosek: Zatem tak jak zauwazono przy histogramach, zadna zmienna nie ma rozkladu normalnego, a przy dopuszczalnym poziomie bledu 5% Test Shapiro-Wilka to udowodnil.

apply(krab[,4:8],2,function(x) silverman.test(x,k=1))

## $Height
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0 
## 
## $Weight
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.09009009 
## 
## $Shucked.Weight
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.5515516 
## 
## $Viscera.Weight
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.1471471 
## 
## $Shell.Weight
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.5955956

Wniosek: Przy dopuszczalnym poziomie bledu 5% test Silvermana pokazal ze rozklad zmiennej Height jest naprawdopodobniej wielomodalny. Pozostale rozklady wskazuja bardziej na teze, ze sa one jednomodalne. Mozna zatem ocenic czy rozklad zmiennej Height jest wiecej niz bimodalny.

nr.modes(hist(krab$Height,plot=F)$counts)

## [1] 11

Wniosek: Oznacza to zatem, ze nie da sie sensownie oszacowac tego parametru, poniewaz rozklad jest za malo zmienny. Zatem zadna zmienna nie jest wielomodalna.

Sprawdzmy zatem ktore zmienne maja obserwacje odstajace. Jest to wazny aspekt przy klastrowaniu.

apply(krab[,-1],2,function(x) grubbs.test(x))

## $Length
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 3.7406, U = 0.9964, p-value = 0.3528
## alternative hypothesis: lowest value 0.1875 is an outlier
## 
## 
## $Diameter
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 3.55873, U = 0.99675, p-value = 0.718
## alternative hypothesis: lowest value 0.1375 is an outlier
## 
## 
## $Height
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 23.58273, U = 0.85707, p-value < 2.2e-16
## alternative hypothesis: highest value 2.825 is an outlier
## 
## 
## $Weight
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 4.06979, U = 0.99574, p-value = 0.09
## alternative hypothesis: highest value 80.10151225 is an outlier
## 
## 
## $Shucked.Weight
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 5.09567, U = 0.99333, p-value = 0.0006477
## alternative hypothesis: highest value 42.184056 is an outlier
## 
## 
## $Viscera.Weight
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 5.28620, U = 0.99282, p-value = 0.0002312
## alternative hypothesis: highest value 21.54562 is an outlier
## 
## 
## $Shell.Weight
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 5.50171, U = 0.99222, p-value = 6.902e-05
## alternative hypothesis: highest value 28.4912475 is an outlier
## 
## 
## $Age
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 5.91289, U = 0.99101, p-value = 6.047e-06
## alternative hypothesis: highest value 29 is an outlier

Wniosek: Przyjmujac kryterium bledu 5% zmienne Height, Shucked.Weight, Viscera.Weight, Shell.Weight i Age maja obserwacje odstajace, podczas gdy pozostale zmienne maja zbyt malo wystarczajacych dowodow na to, ze istnieja w ich przypadku wartosci odstajace.

3. Analiza korelacji i redukcja wymiaru metoda PCA i UMAP

Korelacje policzymy z uzyciem dwoch estymatorow: Pearson’a i Kendall’a aby uniknac sytuacji kiedy mamy zaklocenie zwiazane z brakiem rozkladu normalnego.

heatmaply_cor(cor(krab[,-1],method = 'pearson'))

Wniosek: Analiza korelacji pokazuje, ze wszystkie zmienne sa dodatnio zalezne. Zmienna Age ma najmniejsza korelacje wobec innych zmiennych, a najnizsza korelacje ma z zmienna Shucked.Weight (na poziomie okolo 0.4). W naszym zestawie danych slabiej dodatnio skorelowana jest rowniez zmienna Height z innymi zmiennymi (glownie na poziomie 0.7-0.8). Zmienne Shucked.Weight, Viscera.Weight, Weight, Shell.Weight, Diameter i Length sa silnie dodatnio zalezne. Najsilniejsza korelacje widzimy pomiedzy zmienna Diameter i Length (na poziomie prawie 0.99).

heatmaply_cor(cor(krab[,-1],method = 'kendall'))

Wniosek: Metoda Kendall’a daje bardzo podobne rezultaty do metody Pearson’a. Zmienia sie jedynie wielkosc poziomow korelacji, natomiast jest ona zachowana dokladnie w ten sam sposob.

Redukcja wymiaru PCA:

prcomp(krab[,-1])->pca.krab
summary(pca.krab)

## Importance of components:
##                            PC1     PC2     PC3    PC4     PC5     PC6     PC7
## Standard deviation     16.0148 2.92017 1.34860 0.9263 0.60666 0.13769 0.05570
## Proportion of Variance  0.9567 0.03181 0.00678 0.0032 0.00137 0.00007 0.00001
## Cumulative Proportion   0.9567 0.98856 0.99534 0.9985 0.99991 0.99998 1.00000
##                            PC8
## Standard deviation     0.03057
## Proportion of Variance 0.00000
## Cumulative Proportion  1.00000

df.pca=data.frame(pc1=pca.krab$x[,1],pc2=pca.krab$x[,2],pc3=pca.krab$x[,3],kl=as.factor(krab$Sex))
plot_ly(df.pca,x=~pc1,y=~pc2,z=~pc3,color = ~kl,type='scatter3d')

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Wniosek: Widzimy, ze PC1 ma znaczna wariancje w porownaniu z innymi wartosciami, co oznacza, że zawiera ona najwięcej zmienności danych. PC1 opisuje zdecydowanie wiekszosc danych (96%), a nastepnie razem z PC2 i PC3 jest to zobrazowanie danych na poziomie 99,53%. Na 3 wymiarowej reprezentacji PCA widzimy, ze zadna plec nie jest znaczaco wyrozniona i wyodrebniona od reszty. Mozemy zaobserwowac jedynie ze plec nieokreslona (I) ma niewielkie wyodrebnienie swojej czesci od pozostalych z boku drugiej skladowej glownej.

df.rot<-data.frame(r1=pca.krab$rotation[,1],r2=pca.krab$rotation[,2],r3=pca.krab$rotation[,3])
plot_ly(df.rot,x=~r1,y=~r2,z=~r3)

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

pca.krab$rotation[,1]

##         Length       Diameter         Height         Weight Shucked.Weight 
##    -0.01740306    -0.01438497    -0.00534974    -0.86705518    -0.38119994 
## Viscera.Weight   Shell.Weight            Age 
##    -0.18716833    -0.23523732    -0.10952134

Wniosek: PC1 jest glownie w kierunku przeciwnym do zmiennej Weight, czyli wplywa ona na PC1 w bardzo duzym stopniu. Natomiast zmienna Height ma bardzo niski wplyw na PC1.

Redukcja wymiaru metoda UMAP:

umap.krab<-umap(krab[,-1],n_components = 3)
summary(umap.krab$layout)

##        V1                V2                 V3        
##  Min.   :-9.1023   Min.   :-6.29367   Min.   :-4.529  
##  1st Qu.:-2.9695   1st Qu.:-3.86095   1st Qu.:-3.343  
##  Median : 0.7417   Median :-0.09066   Median :-1.590  
##  Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.000  
##  3rd Qu.: 3.1934   3rd Qu.: 3.67061   3rd Qu.: 3.098  
##  Max.   : 6.0006   Max.   :10.80851   Max.   : 7.511

df.umap=data.frame(um1=umap.krab$layout[,1],um2=umap.krab$layout[,2],um3=umap.krab$layout[,3],kl=as.factor(krab$Sex))
plot_ly(df.umap,x=~um1,y=~um2,z=~um3,color = ~kl,type='scatter3d')

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Wniosek: Po redukcji do trzech wymiarow mozemy zaobserwowac w podsumowaniu jakie wartosci przyjmuje dany wymiar. Wszystkie cechy maja wartosc oczekiwana rowna 0, co oznacza ze wymiar zostal dobrze znormalizowany. Sa one w miare rowno porozdzielane. W porownaniu z PCA, na 3 wymiarowym wykresie UMAP, widzimy juz wieksze wyodrebnienie sie plci nieokreslonej od pozostalych plci krabow. Zatem wskazuje nam to, ze plec nieokreslona (I) bedzie bardziej rozroznialna. Natomiast plcie zenska i meska sa bardzo mocno ze soba pomieszane na wykresie, co wskazuje na mala rozroznialnosc pomiedzy nimi.

4. Analiza skupien - klastrowanie

Aby zaczac klastrowanie warto policzyc entropie:

hist.krab<-lapply(krab[,-1],function(x) hist(x,plot=F))
e.krab<-sapply(hist.krab,function(x) entropy(x$counts))
e.krab

##         Length       Diameter         Height         Weight Shucked.Weight 
##      1.7950410      2.2826319      0.8651742      2.3687285      1.5714585 
## Viscera.Weight   Shell.Weight            Age 
##      1.7894783      2.0379834      1.8319085

Zeby wyznaczyc dobrze entropie naszej zmiennej, na podstawie ktorej mozemy przeprowadzic klastrowanie, trzeba zmienic ja na zmienna o wartosciach liczbowych (klasy 1,2,3).

plec_klasy<-ifelse(krab$Sex == "M", 3,
                   ifelse(krab$Sex == "F", 2,
                          ifelse(krab$Sex == "I", 1, NA)))
entropy(plec_klasy)

## [1] 8.181022

Wniosek: Zmienna Sex (po zamianie wartosci na liczby calkowite) ma wieksza entropie niz pozostale zmienne w zbiorze danych krab, a wiec do klayfikacji trzeba bedzie uzyc algorytmu, ktory wykorzystuje wiele zmiennych.

Zaczniemy od algorytmu k-srednich przy k=3 (bo wiemy ze sa trzy podzialy plci kraba, dzieki czemu mozemy wyznaczyc dany podzial na klasy):

kmeans(krab[,-1], centers=3)->km.krab.3
df.km<-data.frame(x=pca.krab$x[,1],y=pca.krab$x[,2],z=pca.krab$x[,3],type=as.factor(km.krab.3$cluster))
plot_ly(df.km,x=~x,y=~y, z=~z,color=~type,type ='scatter3d')

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Wniosek: Klastrowanie metoda k-srednich jest wzdluz wartosci PC1. Wszystkie trzy klastry sa bardzo blisko siebie, wrecz na siebie zachodza, jednak widzimy dobrze ich wyodrebnienie. Sprawdzmy czy to klastrowanie jest zgodne z podzialem na plec krabow.

table(krab$Sex,km.krab.3$cluster)

##    
##       1   2   3
##   F 373 234 618
##   I  14 972 247
##   M 379 338 718

Wniosek: Do pierwszego klastra nalezy glownie plec nieokreslona (I), jednak nadal mamy wartosci plci zenskiej (F) i meskiej (M). W 2 klastrze mamy juz niejednoznaczny podzial na plec meska i zenska w bardzo podobnych wartosciach, za to malo wartosci plci nieokreslonej. 3 klaster ma lekka przewage plci meskiej, jednak nadal ma duze rozbicie takze na plec zenska. Klastrowanie slabo pokrywa sie z podzialem na plec u krabow, poniewaz wyszlo niejednoznacznie.

Mozemy takze sprobowac metody klastrowania dla duzej ilosci danych - CLARA.

clara(krab[,-1],k=3)->clara.krab.3
df.clara=data.frame(x=pca.krab$x[,1],y=pca.krab$x[,2],z=pca.krab$x[,3],kl=as.factor(clara.krab.3$cluster))
plot_ly(df.clara,x=~x,y=~y, z=~z,color = ~kl,type='scatter3d')

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Wniosek: Klastrowanie metoda clara rowniez dalo nam podzial, w ktorym 3 klastry sa bardzo blisko siebie. Przewage odznacza klasa 3, zas klasa 2 zawiera najmniej wartosci.

Spojrzmy, czy moze to klastrowanie jest zgodne z podzialem na plec krabow:

table(krab$Sex,clara.krab.3$cluster)

##    
##       1   2   3
##   F 451 101 673
##   I 450 727  56
##   M 540 183 712

Wniosek: Uzyskany wynik nadal jest niejednoznaczny i klastrowanie slabo pokrywa sie z podzialem na plec u krabow. W 1 klastrze sa wartosci mocno usrednione, kazda wartosc pasuje prawie tak samo, chociaz plec meska wykazuje przewage. Natomiast w 2 klastrze dopiero mamy dominacje plci nieokreslonej. 3 klaster ma podzial na plec meska i zenska.

Wybor ilosci klastrow (metoda lokcia):

wss<-rep(NA,9)
for(i in 1:9) wss[i]<-kmeans(krab[,-1],centers=i+1)$tot.withinss
plot(1:9,wss,main='Metoda lokcia')

Wniosek: Zatem oznacza to, ze optymalnie jest podzielic zbior danych na 3 lub 4 klastry, poniewaz w tym miejscu krzywa zaczyna wolniej malec.

Mozemy takze zastosowac metode silhouette, aby wybrac ilosc klastrow:

fviz_nbclust(krab[,-1], kmeans, method='silhouette')

Wniosek: Patrzac tutaj na inna metode widzimy, ze jest ona bardziej sklonna ku temu, aby podzielic zbior na 2 klastry.

Klastrowanie hierarchiczne:

plot(hclust(dist(krab[,-1],method = 'minkowski',p=1)))

Wniosek: Typowe klastrowanie hierarchiczne przy tak duzym zbiorze danych nie bedzie optymalnym rozwiazaniem, ze wzgledu na zlozonosc obliczeniowa. Natomiast gdy wykorzystamy bardziej dostosowana i szybsza metode do naszego zbioru danych (hclust), uzyskujemy oczekiwany dendogram. Drzewo hierarchiczne na samym poczatku zostalo podzielone na 2, a pozniej na 3. Nastepnie zaczynaja sie tworzyc zbyt male podgrupy. Oznacza to, ze najlepiej jest wybrac podzial na 3 klastry wedlug tej metody.

5. Klasyfikacja pod nadzorem z uzyciem lasow losowych

rfcv(krab[,-1],as.factor(krab$Sex))->rfcv.krab
rfcv.krab$error.cv

##         8         4         2         1 
## 0.4569741 0.4829181 0.4677627 0.4934498

Wniosek: Mozliwie niski blad walidacji krzyzowej i najbardziej optymalny mamy przy 2 zmiennych. Wszelkie inne podzialy przyniosa nam wiecej bledow. Jednak i tak przy kazdej zmiennej musimy sie liczyc z tym, ze bedzie ona malo miarodajna. Wybierzmy te 2 zmienne.

randomForest(krab[,-1],as.factor(krab$Sex),importance = T,proximity=TRUE)->rf.krab.imp
varImpPlot(rf.krab.imp)

Wniosek: Zmienne sa calkowicie podzielone, jednak kierujac sie Mean Decrease Gini, zmienne Viscera.Weight i Weight sa najlepsze do tego lasu losowego.

randomForest(krab[,c(5,7)],as.factor(krab$Sex))->rf.krab.final
rf.krab.final$confusion

##     F   I   M class.error
## F 474 209 542   0.6130612
## I 161 835 237   0.3227899
## M 509 304 622   0.5665505

Wniosek: Optymalny model lasu losowego zawierajacy 2 zmienne daje maksymalny blad klasyfikacji na poziomie 61%! Oznacza to, ze model lasu losowego ma bardzo niska skutecznosc w tym zestawie danych.

Mozemy takze przeprowadzic analize lasu losowego przy uzyciu danych treningowych i testowych, aby ocenic czy jest to lepszy model, dla ktorego las losowy mialby wyzszy poziom skutecznosci w klasyfikacji danych.

set.seed(12345)
ind <- sample(c(TRUE, FALSE), nrow(krab), replace = TRUE, prob = c(0.7, 0.3))
train <- krab[ind,]
test <- krab[!ind,]

train$Sex<-as.factor(train$Sex)
test$Sex<-as.factor(test$Sex)

rf <- randomForest(Sex~., data=train, proximity=TRUE,mtry=5)
train.tree.pr <- predict(rf, train)
confusionMatrix(train.tree.pr, train$Sex)

## Registered S3 methods overwritten by 'proxy':
##   method               from    
##   print.registry_field registry
##   print.registry_entry registry

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    F    I    M
##          F  837    0    0
##          I    0  860    0
##          M    0    0 1008
## 
## Overall Statistics
##                                      
##                Accuracy : 1          
##                  95% CI : (0.9986, 1)
##     No Information Rate : 0.3726     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 1          
##                                      
##  Mcnemar's Test P-Value : NA         
## 
## Statistics by Class:
## 
##                      Class: F Class: I Class: M
## Sensitivity            1.0000   1.0000   1.0000
## Specificity            1.0000   1.0000   1.0000
## Pos Pred Value         1.0000   1.0000   1.0000
## Neg Pred Value         1.0000   1.0000   1.0000
## Prevalence             0.3094   0.3179   0.3726
## Detection Rate         0.3094   0.3179   0.3726
## Detection Prevalence   0.3094   0.3179   0.3726
## Balanced Accuracy      1.0000   1.0000   1.0000

Wniosek: Las losowy jest przygotowany w dobry sposob na podstawie danych treningowych. Mozemy wiec przejsc do testowania na naszych danych testowych.

test.tree.pr <- predict(rf, test)
confusionMatrix(test.tree.pr, test$Sex)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   F   I   M
##          F 156  22 146
##          I  45 280  69
##          M 187  71 212
## 
## Overall Statistics
##                                          
##                Accuracy : 0.5455         
##                  95% CI : (0.5166, 0.574)
##     No Information Rate : 0.3594         
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.3161         
##                                          
##  Mcnemar's Test P-Value : 0.004697       
## 
## Statistics by Class:
## 
##                      Class: F Class: I Class: M
## Sensitivity            0.4021   0.7507   0.4965
## Specificity            0.7900   0.8601   0.6610
## Pos Pred Value         0.4815   0.7107   0.4511
## Neg Pred Value         0.7315   0.8829   0.7006
## Prevalence             0.3266   0.3140   0.3594
## Detection Rate         0.1313   0.2357   0.1785
## Detection Prevalence   0.2727   0.3316   0.3956
## Balanced Accuracy      0.5960   0.8054   0.5787

Wniosek: Przy uzyciu lasu losowego z wykorzystaniem danych treningowych i testowych uzyskujemy wieksza skutecznosc modelu na poziomie okolo 55%. Najwiecej problemow model wykazuje przy rozroznianiu plci meskiej od zenskiej. Procent skutecznosci w rozroznianiu plci nieokreslonej od pozostalych to 80%.

6. Analiza z uzyciem LDA i QDA

set.seed(12345)
train.ind<-sample(c(TRUE,FALSE),nrow(krab),prob=c(0.7,0.3),replace=T)
train.krab<-krab[train.ind,]
test.krab<-krab[!train.ind,]

lda(Sex~.,data=train.krab)->lda.krab
lda.krab

## Call:
## lda(Sex ~ ., data = train.krab)
## 
## Prior probabilities of groups:
##         F         I         M 
## 0.3094270 0.3179298 0.3726433 
## 
## Group means:
##     Length  Diameter    Height   Weight Shucked.Weight Viscera.Weight
## F 1.446461 1.1359319 0.3927569 29.63187      12.626823       6.506329
## I 1.075828 0.8190116 0.2720349 12.44260       5.503116       2.655310
## M 1.396987 1.0915303 0.3757813 27.68066      12.016771       6.035406
##   Shell.Weight      Age
## F     8.546578 11.14337
## I     3.721927  7.87907
## M     7.895898 10.75397
## 
## Coefficients of linear discriminants:
##                         LD1        LD2
## Length          3.970582757  3.1115703
## Diameter       -5.552422675 -8.8736382
## Height         -3.262518773 -2.5541100
## Weight         -0.014131327 -0.1134331
## Shucked.Weight -0.009340328  0.5710888
## Viscera.Weight -0.197865424 -0.1269375
## Shell.Weight    0.096833756 -0.1154048
## Age            -0.127632669  0.2478620
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9883 0.0117

train.krab$Sex<-as.factor(train.krab$Sex)
plot(lda.krab, col=as.numeric(train.krab$Sex))

Wniosek: Prawdopodobienstwa wystapienia danej grupy plci sa na podobnym poziomie. Mozemy takze wyczytac, ze dla plci nieokreslonej (I) mamy najmniejsze wartosci w podsumowaniu ich srednich. Wplyw na klasyfikacje ma przede wszystkim LD1 z cecha Length. Na podstawie wykresu widzimy ze zmienne mocno zachodza na siebie i nie sa dobrze rozroznialne.

Popatrzmy na jakosc klasyfikacji na wybranych parach zmiennych:

partimat(Sex ~ Length + Diameter + Height + Age, data=train.krab, method="lda")

partimat(Sex ~ Weight + Shucked.Weight + Viscera.Weight + Shell.Weight, data=train.krab,method="lda")

Wniosek: Jakosc klasyfikacji jest na dosc niskim poziomie, poniewaz nawet dla zmiennych o duzej korelacji wykazuje on poziom bledu na poziomie okolo 48%.

Sprawdzmy teraz calkowite sprawdzenie jakosci klasyfikacji z uzyciem LDA:

train.krab$Sex<-as.factor(train.krab$Sex)
predict(lda.krab,train.krab)->pre.train
cat("Dane treningowe:","\n")

## Dane treningowe:

table(train.krab$Sex,pre.train$class)

##    
##       F   I   M
##   F 258 111 468
##   I  30 657 173
##   M 209 208 591

predict(lda.krab,test.krab)->pre.test
cat("Dane testowe:","\n")

## Dane testowe:

table(test.krab$Sex,pre.test$class)

##    
##       F   I   M
##   F 122  53 213
##   I   5 287  81
##   M 102  77 248

Wniosek: Model w niewielkiej wiekszosci przypadkow poprawnie ocenil plec krabow, jednak czesto mocno mylil sie pomiedzy plcia meska (M) i zenska (F). Najlepiej wyszla klasyfikacja dla plci nieokreslonej (I).

Sprawdzmy teraz analize z uzyciem QDA:

set.seed(12345)
train.ind.q<-sample(c(TRUE,FALSE),nrow(krab),prob=c(0.7,0.3),replace=T)
train.krab.q<-krab[train.ind.q,]
test.krab.q<-krab[!train.ind.q,]
qda(Sex~.,data=train.krab.q)->qda.krab
qda.krab

## Call:
## qda(Sex ~ ., data = train.krab.q)
## 
## Prior probabilities of groups:
##         F         I         M 
## 0.3094270 0.3179298 0.3726433 
## 
## Group means:
##     Length  Diameter    Height   Weight Shucked.Weight Viscera.Weight
## F 1.446461 1.1359319 0.3927569 29.63187      12.626823       6.506329
## I 1.075828 0.8190116 0.2720349 12.44260       5.503116       2.655310
## M 1.396987 1.0915303 0.3757813 27.68066      12.016771       6.035406
##   Shell.Weight      Age
## F     8.546578 11.14337
## I     3.721927  7.87907
## M     7.895898 10.75397

Wniosek: Rozlozenie prawdopodobienstw jest takie samo jak w przypadku LDA.

Sprawdzmy macierz bledu dla QDA:

train.krab.q$Sex<-as.factor(train.krab.q$Sex)
predict(qda.krab,train.krab.q)->pre.train.q
cat("Dane treningowe:","\n")

## Dane treningowe:

table(train.krab.q$Sex,pre.train.q$class)

##    
##       F   I   M
##   F 384 196 257
##   I  62 731  67
##   M 308 323 377

predict(qda.krab,test.krab.q)->pre.test.q
cat("Dane testowe:","\n")

## Dane testowe:

table(test.krab.q$Sex,pre.test.q$class)

##    
##       F   I   M
##   F 168  97 123
##   I  25 316  32
##   M 154 123 150

Wniosek: W przypadku QDA mamy troche inaczej rozlozone nasze wyniki danych testowych. Mianowicie o wiele wiecej sklasyfikowano plci nieokreslonej w poprawny sposob, jednak nadal widzimy mocne podobienstwa miedzy klasyfikacja plci meskiej i zenskiej. Sa minimalne przewagi w poprawnej klasyfikacji, jednak nadal jest ona na slabym poziomie.

Wnioski koncowe

Na podstawie danych zgromadzonych na temat krabow hodowlanych i ich cech fizycznych moglismy stwierdzic wiele zaleznosci. Przede wszystkim zmienna Weight wyroznila sie znaczaca zmiennoscia oraz duzym wplywem na analizowane dane. Zmienne Diameter oraz Length sa ze soba najbardziej skorelowane i najbardziej na siebie oddzialuja. Ale warto zaznaczyc, ze tak naprawde wszystkie zmienne posiadaja znaczna korelacje i na siebie w jakims stopniu oddzialuja. Najmniejsza korelacje ma zmienna Age i oznacza to, ze wiek u krabow nie jest az tak zalezny od innych czynnikow. Za to, wszystkie zmienne dotyczace wagi sa miedzy soba mocno skorelowane. Przy zmiennej Height wiemy za to, ze posiada ona pewne wartosci odstajace, ktore wplywaja na jej asymetrie i na tworzenie roznych modeli klasyfikacyjnych, poniewaz w jakims stopniu zakloca ona wyniki. Rozklady danych sa jednomodalne, zatem wiekszosc danych koncentruje sie wokol jednej wartosci.

Przy probach redukcji wymiaru metoda PCA oraz UMAP widoczne jest jak bardzo plec meska i zenska u krabow nie jest az tak bardzo rozroznialna i zadna cecha ich miedzy soba nie rozroznia. Natomiast moglismy zauwazyc na podstawie kolejnych analiz, ze to plec nieokreslona byla najlepsza do wytypowania przez modele, poniewaz mozna ja bylo najlatwiej rozroznic. Na podstawie danych prob klastrowania widzimy, ze sie one od siebie znacznie roznia i kazda metoda dala rozny wniosek. Oznacza to zatem, ze optymalny numer klastrow, w tym modelu, waha sie pomiedzy 2 a 4. Przy zastosowaniu metody klasyfikacji pod nadzorem z wykorzystaniem lasu losowego uzyskalismy niepowodzenie, poniewaz wynikiem koncowym bylo 61% bledu klasyfikacyjnego. Natomiast przy doborze danych treningowych oraz testowych dalo sie zoptymalizowac ten blad do 45%. Oznacza to zatem, ze ten typ modelu nie jest dobry do tego zbioru danych. Czesto wystepowaly problemy z klasyfikacja plci meskiej i zenskiej w poprawny sposob.

Przy uzyciu metody analizy z wykorzystaniem LDA i QDA dochodzimy do tych samych wnioskow, ktore mowia nam o tym, ze model czesto myli sie pomiedzy plciami zenska i meska, jednak rozroznial plec nieokreslona na o wiele wyzszym poziomie. Zatem kraby o plci meskiej i zenskiej sie od siebie, az tak nie roznia, natomiast gdy dochodzi do pomiarow krabow o plci nieokreslonej to juz latwiej nam rozroznic o jakie wyniki moze chodzic.

Rezultaty te sugeruja potrzebe dalszych badan w celu optymalizacji modeli klasyfikacyjnych oraz poszukiwania dodatkowych cech, ktore moglyby lepiej rozrozniac plec krabow hodowlanych.